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We study the behavior of cold nuclear matter near saturation density (po) and very low tem- 
perature using classical molecular dynamics. We used three different (classical) nuclear interaction 
models that yield 'medium' or 'stiff' compressibilities. For high densities and for every model the 
ground state is a classical crystalline solid, but each one with a different structure. At subsaturation 
densities, we found that for every model the transition from uniform (crystal) to non-uniform matter 
occurs at p « 0.12/ra -3 = 0.75/>o- Surprisingly, at the non-uniform phase, the three models produce 
'pasta-like' structures as those allegedly present in neutron star matter but without the long-range 
Coulomb interaction and with different length scales. 
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I. INTRODUCTION 

The study of cold nuclear matter at subsaturation den- 
sities is important for a variety of topics which include the 
structure of neutron star crusts, the equation of state of 
nuclear matter, as well as many other phenomena related 
to the study of heavy ion reactions. Pioneering stud- 
ies based on the compressible liquid drop model p], Q, 
Hartree-Fock method Q , and energy minimization tech- 
niques 0, [H showed that transitions from spherical nuclei 
to shapes such as rods, slabs, tubes, spherical bubbles, 
are to be expected in cold symmetric nuclear matter em- 
bedded in an electron gas. More recent investigations 0- 
have used dynamical methods to study such transi- 
tions in neutron star crust environments. 

A major component in all of these studies has been the 
presence of an all embedding electron gas which, in neu- 
tron stars, is expected to be produced by an abundance 
of /3-decays. In fact, a Coulomb interaction has been 
thought of as an essential ingredient for the formation of 
the rich "pasta" like structures. Koonin [5], for instance, 
explains the transitions between different topologies in 
terms of a competition between a short-range nuclear 
surface energy which becomes minimized through aggre- 
gation, and a long-range Coulomb energy which gets re- 
duced by an opposite dispersion. More recent studies of 
Horowitz and coworkers [16j re-examine this assertion as 
an example of frustration, a phenomenon that emerges 
from the impossibility to simultaneously minimize all in- 
teractions, and yields a large number of low-energy con- 
figurations. Indeed in all of the previous studies listed, 
the Coulomb interactions have been approximated ei- 
ther by electrically grounded surfaces Q , uniform electric 
charge densities p , Thomas- Fermi screened Coulomb po- 
tential (see, e.g. or by an Ewald summation [l8| . 

The importance the Coulomb repulsion has on the for- 
mation of the pasta-like structures prompted our previ- 
ous study [l9j] which dissected the role of this interac- 



tion as a function of density, temperature and isotopic 
content through the use of classical molecular dynamics 
simulations. In an unexpected outcome, however, such 
study, which varied the strength of the electric interac- 
tion from full to none, observed that the pasta structures, 
namely, "gnocchi" , "spaghetti" , "lasagna" and their anti- 
structures (i.e. those obtained by replacing particles by 
holes and viceversa) seemed to exist even in the absence 
of any Coulomb interaction. Although the structures ap- 
peared modified somewhat in their overall scales, topol- 
ogy and location in the density-temperature plane, the 
attractive and repulsive components of the nuclear po- 
tential were found to be sufficient to give rise to a rich 
pasta-like structure in infinite nuclear systems at subsat- 
uration densities and low temperatures. 

The implications of these findings on the nuclear equa- 
tion of state are immediate. The fact that nuclear matter 
at subsaturation densities (but non-zero temperatures, 
T > 1 MeV) decomposes into a mix of liquid and gaseous 
phases [20[, while at lower temperatures (T < 1 MeV) 
it self- assembles into pasta- like objects, gives rise to very 
interesting questions that we would like to address in the 
present work. 

Explicitly, our interest is in understanding the struc- 
ture of infinite cold nuclear matter at subsaturation den- 
sities. At what density does an infinite cold nuclear sys- 
tem abandons a uniform phase (indeed crystalline at zero 
temperature) to produce pasta-like objects? What struc- 
tures are formed at various densities? Other related ques- 
tions are the dependence of this phenomena on the par- 
ticular inter-particle potentials used. 

The approach to be followed is a combination of classi- 
cal molecular dynamics, fragment recognition algorithms 
and topological analysis tools to study the structure of 
infinite nuclear matter at subsaturation densities and 
very cold temperatures. In the next section the clas- 
sical molecular dynamics model used will be briefly re- 
viewed for completeness. Section [Till presents a detailed 
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study of the structure of symmetric nuclear matter with 
medium compressibility, followed in Section IIVI by simi- 
lar studies with other potentials that have been used in 
the past to study nuclear systems. We present a sum- 
mary of our main results in Section IVl and outline our fu- 
ture study that will explore the connection between our 
findings with those at higher temperatures, i.e. the link- 
age from the nuclear pasta to the nuclear fragmentation 
regime. 



II. CLASSICAL MOLECULAR DYNAMICS 

We use a classical molecular dynamics (CMD) model 
to study infinite nuclear matter at low temperatures and 
subsaturation densities; the use of molecular dynamics 
to study nuclear reactions was pioneered by Wilets and 
coworkers (2l| and advanced by Pandharipande [22| and 
others [H, Recently, classical molecular dynamics 
models have been used to study cold nuclear matter in 
neutron star crusts environments [H, l25l-[27| ; in particu- 
lar the CMD model, which was developed to study nu- 
clear reactions |28l-l37| , has been adapted to study infinite 
nuclear systems under such conditions [HI, [l9| . 

In our study, the trajectories of the nucleons are gov- 
erned by classical equations of motion dictated by forces 
produced by the Pandharipande [22| potentials: 

V np (r) = V r [exp(-p r r)/r - exp(-p r r c )/r c ] 

- V a [exp(-p a r)/r - exp(-p a r c )/r c ] 
V NN (r) = V [exp(-fjL r)/r - exp(-por c )/r c } , 

where the attractive potential between a neutron and a 
proton is V npi and the repulsive interaction between nn 
or pp is Vat at; they both use a cutoff radius of r c = 5.4 
fm after which the potentials are set to zero. The 
Yukawa parameters /i r , \i a and po were phenomenolog- 
ically adjusted by Pandharipande to yield a saturation 
density of po = 0.16 fm~ 3 , a binding energy E(po) = 
16 MeV/nucleon and a compressibility listed in [22j as 
250 MeV for the "Medium" model, and 535 MeV for 
the "Stiff"' 0. 

The trajectories of all nucleons are obtained solving 
the classical equations of motion using a Verlet algorithm 
with appropriate energy conservation. To mimic an in- 
finite system thousands of nucleons were placed in cubic 
cells under periodic boundary conditions. The densities 
were enforced by placing either 1000 protons and 1000 
neutrons (to study systems symmetric in isospin, x = 
z/A = 0.5) or 1000 protons and 2000 neutrons (for isospin 
asymmetric cases, x = 0.3) in cubical boxes with sizes ad- 
justed to have densities between 0.001 fm~ 3 < p < po; 
previous studies already presented samples of structures 
obtained through this method with [l5| and without [l9| 
an embedding electron gas. To study systems at differ- 
ent temperatures, the nuclear matter is force- heated or 
cooled using the Andersen thermostat procedure [38| to 
control the temperature. 




FIG. 1. Binding energies per nucleon for systems obtained 
with the Pandharipande medium potential for crystalline lat- 
tices with simple cubic (SC), face-centered cubic (FCC) and 
body-centered cubic (BCC) geometries, and using molecular 
dynamics at T = 0.001 MeV (CMD). The structures cor- 
responding to the four labeled points ("A" through "D") are 
shown in Figure [2] 



CMD provides the position and momentum of all 
nucleons at all times which, when coupled to a clus- 
ter recognition algorithm allows the characterization of 
the nuclear structure. In the case of very low tempera- 
tures, the identification of clusters can be achieved with 
a "Minimum Spanning Tree" (MST) algorithm [39|,|40|; 
which identifies clusters simply by finding nucleons that 
are closer to each other than some clusterization radius, 
which we set to 3.0 fm. The MST method used was 
adapted to recognize fragments that extend into any of 
the adjacent periodic cells. 

To simulate cold nuclear matter in the liquid-gas co- 
existence region we operate in the range of densities of 
0.01 fm~ 3 < p < p . To examine the behavior at very 
cold temperatures we focus on < T < 1.0 MeV . At 
a difference from our previous studies [HI, [HI in which 
configurations were sampled hundreds of times in an sta- 
tionary ergodic process, at these low temperatures it is 
not necessary to perform such sampling as the produced 
configurations are either crystalline or have very low mo- 
bilities. 

Our procedure is two fold. To study the uniform phase 
exactly at zero temperature, a crystalline structure at a 
given density is adopted and its energy per nucleon cal- 
culated by direct summation between all nucleons. The 
structure then is scaled out to reach a different density 
where the energy is evaluated again, and so on. This pro- 
cedure, identical to the one used by Pandharipande [22| 
and Koonin [E| produces the characteristic "U" shapes 
when the binding energy is plotted versus the density, 
with a minimum signalling the saturation, or normal, 
nuclear density. Figure [1] shows such curve for the Pand- 
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FIG. 2. (Color online) Structures corresponding to the la- 
beled points of Figure [1] Point A corresponds to a formation 
in the regular lattice, while the rest of the points are pasta 
structures. 



FIG. 4. (Color online) Structures corresponding to the labeled 
points of Figure [3] obtained with the Pandharipande medium 
potential at T = 1.0 MeV. 



haripande medium potential. 

The second method uses CMD and starts from a ran- 
dom positioning of the nucleons in the central simple cu- 
bic cell (all replicated in the surrounding 26 cells) and en- 
dowed with some velocity distribution corresponding to a 
given initial temperature. The system is then heated up 
to a large temperature (T > 1 MeV) and then brought 
down to the final desired temperature using the Ander- 
sen thermostat procedure in small temperature steps. We 
take T = 0.001 MeV as zero temperature, other values 
explored are in the range 0.001 MeV < T < 1. After 
reaching equilibrium, the clusters are identified and the 
analysis tools described in [15] are used to visualize and 
characterize the produced structures. It is with these 
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FIG. 3. Binding energies per nucleon for systems obtained 
with the Pandharipande medium potential at the listed tem- 
peratures. 



tools that we now proceed to study cold nuclear mat- 
ter under a variety of conditions and for several types of 
potentials. 



III. MEDIUM COMPRESSIBILITY MATTER 

Figure [T] shows the CMD zero temperature energies 
for the symmetric x = 0.5 case with the Pandharipande 
medium potential. As can be seen, CMD agrees with in- 
superable precision the simple cubic lattice calculations 
up to a density of p ~ 0.13 /m -3 , for lower densities the 
systems breaks -not in a liquid-gas phase, as expected 
for higher temperatures, but- into pasta-like objects of 
different shapes. Also shown in figure [1] are the ener- 
gies expected for a uniform face centered cubic structure 
which, being higher in energy, do not correspond to the 
T = case. 

As an aside, the SC curve of Figure [1] can be used 
to fit a polynomial fit around its minimum to extract 
an analytic expression which can then be used to calcu- 
late the compressibility at around the saturation density 
through K = 9pl[d 2 (E/A)/dp 2 ] Po . In the case of the 
medium Pandharipande potential we find its value to be 
197 MeV, much less than the value of 250 MeV listed 
by its creators. 

Figure [2] shows the structures corresponding to the four 
densities labeled from "A" to "D" in Figure [TJ These re- 
sults are not entirely new, in fact they coincide with those 
obtained by Williams and Koonin in 1985 with a static 
procedure J5[ . They began by describing the matter as a 
periodic simple cubic lattice at a given number density of 
nucleons, the system was then allowed relax freely until 
a local minimum was achieved; although the method did 
not take into account the possibility of having free nucle- 
ons and used a fixed geometry by construction, it yielded 
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FIG. 5. Curvature - Euler coordinates of the structures of 
Figure [3] The lines connect points with the same densities 
but temperatures varying from T — 0.001 MeV to 1.0 MeV . 



TABLE I. Classification Curvature - Euler 
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Anti-Gnocchi 
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results comparable to those obtained by our dynamical 
method. For instance, structure C in Figure [2] -which 
corresponds to what in [l5[ was dubbed as "lasagnas"- 
was also found by [E[ as "alternating slabs of matter and 
vacuum" . Structure D is also very close to Koonin's mix 
of three-dimensional nuclei. It was a case of good judge- 
ment that Williams and Koonin chose a simple cubic as 
their basis, had they preferred to use FCC or BCC in- 
stead, their saturation density and binding energy would 
have been off from the correct values. Likewise, their as- 
sumption of no free neutrons turned out to be correct for 
the case of x = 0.5. 

Continuing with the study at higher temperatures, Fig- 
ure [3] shows the same type of results as Figure [T] for tem- 
peratures 0.001 MeV < T < 1.0 MeV. As it can be 
seen in the figure, at densities p > 0.12 fm~ 3 and at all 
of these temperatures the curves follow the "U" shape 
characteristic of the uniform T = crystalline phase. As 
in the case of zero temperature, if the density decreases 
below, say p < 0.12 fm~ 3 the systems again move away 
from the uniform phase forming complex arrangements. 
Figure [H shows the structures obtained for the highest 
temperature studied (T = 1.0 MeV) for the four densi- 
ties labeled from "a" to "d" in Figure [3l 

These non- homogeneous structures being formed in the 
low density region can be characterized using the mean 
curvature and Euler number. As explained in detail 
in [l5|, different structures have distinct values of these 
variables and, in general, follow the pattern outlined in 



FIG. 6. Binding energies per nucleon for systems obtained 
with the Simple Semiclassical Potential for crystalline lat- 
tices with simple cubic (SC), face-centered cubic (FCC) and 
body-centered cubic (BCC) geometries, and using molecular 
dynamics at T = 0.001 MeV (CMD). The structures cor- 
responding to the four labeled points ("A" through "D") are 
shown in Figure [3 

Table HJ as a reference, the perfect crystals formed at 
T = and p > 0.12/m -3 (v.g. point "A" in Figures [1] 
and [3]) all have zero mean curvature and zero Euler num- 
ber. 

The curvature-Euler coordinates of the labeled struc- 
tures of Figure [3] are presented in Figure [5j points joined 
by lines all have the same densities but their tempera- 
tures vary through T = 0.001, 0.1, 0.5, 0.6, 0.8 and 
1.0 MeV with the letter labels corresponding to those of 
Figure 03 Again, the cases "A- a" at normal density (po) 
correspond to crystalline structures and all have zero cur- 
vature and Euler number at the temperatures studied. 
The case "B-b" of density 0.01 fm~ 3 goes from being a 
slab with "anti-spaghetti" (tubes of empty space) to a 
slab with "anti-gnocchis" (spherical holes). As pointed 
out before, "C" goes from being a perfect "lasagna" to a 
sloppy slab filled with holes ( "anti-gnocchis" ). Structures 
"D-d" , on the other side, go from spherical "gnocchis" to 
deformed "gnocchis" in what can be considered as the 
beginning of the transition from the pasta phase to the 
liquid-gas mixed phase which is known to exist at higher 
temperatures. 



IV. HIGHER COMPRESSIBILITY 
POTENTIALS 

As the comportment presented in the previous section 
is bound to be potential-dependent, it is instructive to 
repeat the study using other nuclear interactions to ex- 
tract generalities of the behavior of nuclear matter at low 
temperature and subsaturation densities. In particular 
we explore potentials with higher values of the compress- 
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FIG. 7. (Color online) Structures obtained with the SSP 
corresponding to the labeled points of Figure [6] Again, points 
A and B correspond to formations in the regular lattice, while 
points C and D are pasta structures. 

ibility that have been used in the past to study nuclear 
matter. 



A. A simple semiclassical potential 

A higher compressibility potential that has been used 
for a variety of studies of nuclear matter [25l-[27| is the one 
described by its creators [16] as a "simple semiclassical 
potential" (SSP). In summary, the SSP is composed of: 

V np (r)=ae- r2 / A + [b-c]e- r2 / 2A , 

V NN (r) = ae~ r ^ A + [b + c)e~ r ' ^ 2A , 

where, again, the potential between a neutron and a pro- 
ton is attractive, and that between like particles is repul- 
sive. The parameters a, b, c, and A have been adjusted 
to have the proper energy and density scales to mimic 
nuclear matter. Notice that, at a difference from the 
Pandharipande potentials, the SSP do not have repul- 
sive hard cores. As shown in the cited literature, these 
potentials, when coupled with a screened Coulomb inter- 
action, reproduce the expected long range structure of 
the nuclear pasta phases. 

Figure [6] shows the binding energies per nucleon ob- 
tained from symmetric (x = 0.5) systems constructed 
with the SSP in crystalline lattices with simple cubic 
(SC), face-centered cubic (FCC) and body-centered cu- 
bic (BCC) geometries. When used in our molecular 
dynamics code at T = 0.001 MeV -and without the 
screened Coulomb potential-, the SSP produces struc- 
tures with the energies labeled in Figure [6] as CMD. 

As in the case of the medium Pandharipande poten- 
tial, at densities p > 0.12 fm~ 3 the CMD results agree 
with the lowest-energy crystalline structure, which in this 
case is the BCC, and at lower densities the systems form 



pasta-like objects. The structures corresponding to the 
four labeled points ( "A" through "D" ) are shown in Fig- 
ure [3 

Noticeably, the "normal" density point of this poten- 
tial is at the correct value (i.e. at p = 0.16 fm~ 3 ) but 
at a slightly lower binding energy of about —20.3 MeV. 
Again, using a polynomial fit to the bottom part of the 
"U" of the BCC curve of Figure [6j allows us to extract 
a value of the compressibility of the order of 290 MeV, 
much higher than the value we obtained for the Pand- 
haripande medium potential. 4 

Except for minor differences, the scenario that the 
SSP presents is very similar to that obtained with the 
Pandharipande medium potential. Namely, it is possible 
to obtain rich pasta-like structures in an SSP medium 
at cold temperatures and subsaturation densities without 
the modulating effect of the Coulomb interaction. 



B. Stiff Pandharipande potential 

We now repeat the same calculations for the stiff ver- 
sion of the Pandharipande potential; see [22| for the 
proper Yukawa parameters. This interaction produces 
the E/A versus p shown in Figure[8]for symmetric nuclear 
matter and zero temperature in simple cubic (SC), face- 
centered cubic (FCC) and body-centered cubic (BCC) 
crystals. Also presented are the binding energies of the 
structures obtained with CMD at T = 0.001 MeV. 

The authors of the potential chose the prameters to 
correctly reproduce the binding energy at the normal 
density of p = 0.16 fm—3 assuming a SC structure. 
Interestingly, the CMD simulations tell us that in the 
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FIG. 8. Binding energies per nucleon for systems obtained 
with the Pandharipande stiff potential for crystalline lat- 
tices with simple cubic (SC) and face-centered cubic (FCC), 
and using molecular dynamics at T = 0.001 MeV (CMD). 
The structures corresponding to the four labeled points ( "A" 
through "D") are shown in Figure [9] 
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Again, except for minute differences, stiff nuclear 
matter also produces rich pasta-like structures at cold 
temperatures and subsaturation densities without any 
Coulomb interaction. 



FIG. 9. (Color online) Structures corresponding to the labeled 
points of Figure [8] Again, points A and B correspond to 
formations in the regular lattice, while points C and D are 
pasta structures. 



range 0.13 < p < 0.17 fm~ 3 the preferred crystalline 
structure is FCC and it is only at p > 0.17 MeV that 
it becomes SC. 

For this potential, the true saturation point for cold 
matter occurs at p = 0.15 fm~ 3 and with a binding 
energy of about 16.5 MeV. Using the FCC curve to 
estimate the compressibility at the saturation point as 
done with the other two potentials yields a value of about 
366 MeV, jnuch smaller than the value of 535 MeV re- 
ported in (22j . 

The departure from the homogeneous phase into the 
pasta structures is apparent at around p ~ 0.13 fm~ 3 ; 
the structures corresponding to the four labeled points 
( "A" through "D" ) are shown in Figure 



V. CONCLUDING REMARKS 

We studied the behavior of cold symmetric nuclear 
matter in a wide range of densities using three different 
classical models of nuclear interaction. 

The most striking result is the existence of pasta-like 
strutures in absence of the Coulomb interaction. It was 
believed that the long range repulsive Coulomb inter- 
action was needed for the formation of these pasta-like 
structures, but from our results it seems that is not the 
case. 

Apparently, the Coulomb interaction only serves to ex- 
pand the menu of possible structures, and change their 
relative length scales. 

The general picture for all the potentials is the same: 
Above a saturation density the system presents a uniform 
phase, and below it the system forms a variety of pasta- 
like structures. 

The overall shape and types of pasta found are indeed 
model-dependent but at zero temperature the density at 
which the transition from uniform to pasta occurs is the 
same for all potentials. This is likely to happen because 
all the potentials have the same range. 
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